from pyqubo import Array
import itertools
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from haversine import haversine, Unit
plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号
import time
import math
import pyqubo
from openjij import SQASampler
import math

start_time = time.time()  #开始时间
num_points = 111
coordinates = np.random.rand(num_points, 2)  # 生成0到1之间的随机数
# 转换为DataFrame格式，和原来的文件形式保持一致
df = pd.DataFrame(coordinates, columns=[0, 1])

t = 5
g_num = 4#产生地数
c_num = 3#收集中心数
r_num = 2#回收中心数目
fix_c = 1000 #收集中心固定费用
fix_r = 2000 #回收中心固定费用
fix_k = 1000#一级车辆的固定费用
fix_l=1000#二级车辆的固定费用
fix_z=5000#二级车辆的固定费用
T_k=100#一级车辆的运输费用
T_l=100#二级车辆的运输费用
T_z=100
T_El=1000#二级车辆的运输固废的费用
E_k=500
E_l=500
E_z=1500
lenk_c = 2
k_c=lenk_c*c_num
lenl_r=2
l_r=lenl_r*r_num
lenz_r=2
z_r=lenz_r*r_num

GM = [i for i in range (0,g_num)]#固废产生地列表
CM = [i for i in range (g_num,g_num+c_num)]#收集中心列表
RM = [i for i in range (g_num+c_num,g_num+c_num+r_num)]#回收中心列表
T = [i for i in range (t)]
T1 = [i for i in range (t-1)]

N=GM+CM+RM
CGM = CM+GM
CRM = CM+RM
GRM = GM+RM
GCRM = CM+RM+GM

#读取txt文件
x_loc = list(df.iloc[:,1])
y_loc = list(df.iloc[:,0])

# #画出坐标散点图
# plt.scatter(x_loc[0:g_num],y_loc[0:g_num],c='r')
# plt.scatter(x_loc[g_num:g_num+c_num],y_loc[g_num:g_num+c_num],c='b',marker='d')
# plt.scatter(x_loc[g_num+c_num:g_num+c_num+r_num],y_loc[g_num+c_num:g_num+c_num+r_num],c='g',marker='p')
#
# for i in GM:
#     plt.annotate('产生地/需求地''%d'%i,(x_loc[i]+0.001,y_loc[i]))
# for i in CM:
#     plt.annotate('一级中心''%d'%i,(x_loc[i]+0.001,y_loc[i]))
# for i in RM:
#     plt.annotate('二级中心''%d'%i,(x_loc[i]+0.001,y_loc[i]))
# plt.show()

#各个节点之间的弧
A = [(i,j) for i in N for j in N if j!=i]

#各节点之间的距离
C = {(i, j): round(haversine(( y_loc[i],x_loc[i]), (y_loc[j],x_loc[j]), unit=Unit.KILOMETERS), 1) for i in N for j in N if j != i}

#初始化V_ijt 和 Ce_rt，作为决策变量s和x字典的键值
K=[i for i in range(k_c)]
L=[i for i in range(k_c,k_c+l_r)]
Z=[i for i in range(k_c+l_r,k_c+l_r+z_r)]
V=K+L+Z
a = 1+math.ceil(math.log2(lenk_c))
b = 1+math.ceil(math.log2(lenl_r))
c = 1+math.ceil(math.log2(lenz_r))
shape_itm = (len(N), len(T), len(V))

X_i_t_m = Array.create('X_i_t_m', shape=shape_itm, vartype='BINARY')
X_i = Array.create('X_i', shape=(len(CRM)), vartype='BINARY')  # 直接使用CM的长度
V_j = Array.create('V_j', shape=(len(V)), vartype='BINARY')
a1 = Array.create('a1', shape=(a,c_num), vartype='BINARY')
b2 = Array.create('b2', shape=(b,r_num), vartype='BINARY')
c3 = Array.create('c3', shape=(c,r_num), vartype='BINARY')
d4 = Array.create('d4', shape=(1,k_c,t), vartype='BINARY')
e5 = Array.create('e5', shape=(1,l_r,t), vartype='BINARY')
f6 = Array.create('f6', shape=(1,z_r,t), vartype='BINARY')

def build_objective(X_i: pyqubo.array.Array,X_i_t_m: pyqubo.array.Array,V_j: pyqubo.array.Array) :
    H = (
            sum(fix_c * X_i[i - len(GM)] for i in CM)
            + sum(fix_k * V_j[i] for i in K)
            + sum(
        C[i, j] * X_i_t_m[i, t, m] * X_i_t_m[j, t + 1, m] * T_k for i in CGM for j in CGM for t in T1 for m in K if
        j != i if j != i)

            + sum(fix_r * X_i[r - len(GM)] for r in RM)
            + sum(fix_l * V_j[m] for m in L)
            + sum(fix_z * V_j[m] for m in Z)

            + sum(
        C[i, j] * X_i_t_m[i, t, m] * X_i_t_m[j, t + 1, m] * T_l for i in CRM for j in CRM for t in T1 for m in L if
        j != i)
            + sum(
        C[i, j] * X_i_t_m[i, t, m] * X_i_t_m[j, t + 1, m] * T_z for i in GRM for j in GRM for t in T1 for m in Z if
        j != i)

            + sum(
        C[i, j] * X_i_t_m[i, t, m] * X_i_t_m[j, t + 1, m] * E_k for i in CGM for j in CGM for t in T1 for m in K if
        j != i)
            + sum(
        C[i, j] * X_i_t_m[i, t, m] * X_i_t_m[j, t + 1, m] * E_l for i in CRM for j in CRM for t in T1 for m in L if
        j != i)
            + sum(
        C[i, j] * X_i_t_m[i, t, m] * X_i_t_m[j, t + 1, m] * E_z for i in GRM for j in GRM for t in T1 for m in Z if
        j != i)

    )
    return H

def build_rule1(X_i_t_m: pyqubo.array.Array,X_i: pyqubo.array.Array,a1: pyqubo.array.Array) :
  H = pyqubo.Constraint(
      sum(
          (sum(X_i_t_m[i, 0, m] for m in K)-lenk_c*X_i[i-g_num] + sum((2**l)*a1[l,i-g_num] for l in range(a)))**2
          for i in CM),
          'w_1')
  return H

def build_rule2(X_i_t_m: pyqubo.array.Array,X_i: pyqubo.array.Array,b2: pyqubo.array.Array) :
  H = pyqubo.Constraint(
      sum(
          (sum(X_i_t_m[i, 0, m] for m in L)-lenl_r*X_i[i-g_num] + sum((2**l)*b2[l,i-g_num-c_num] for l in range(b)))**2
          for i in RM),
          'w_1')
  return H


def build_rule3(X_i_t_m: pyqubo.array.Array,X_i: pyqubo.array.Array,c3: pyqubo.array.Array) :
  H = pyqubo.Constraint(
      sum(
          (sum(X_i_t_m[i, 0, m] for m in Z)-lenz_r*X_i[i-g_num] + sum(2**l*c3[l,i-g_num-c_num] for l in range(c)))**2
          for i in RM),
          'w_1')
  return H

def build_rule4(X_i_t_m: pyqubo.array.Array,V_j: pyqubo.array.Array) :
  H = pyqubo.Constraint(
      sum(
          (sum(X_i_t_m[i, 0, m] for i in CM) -V_j[m]) **2
          for m in K),
          'w_1')
  return H

def build_rule5(X_i_t_m: pyqubo.array.Array,V_j: pyqubo.array.Array) :
  H = pyqubo.Constraint(
      sum(
          (sum(X_i_t_m[i, 0, m] for i in RM) -V_j[m]) **2
          for m in L),
          'w_1')
  return H

def build_rule6(X_i_t_m: pyqubo.array.Array,V_j: pyqubo.array.Array) :
  H = pyqubo.Constraint(
      sum(
          (sum(X_i_t_m[i, 0, m] for i in RM) -V_j[m]) **2
          for m in Z),
          'w_1')
  return H

def build_rule7(X_i_t_m: pyqubo.array.Array) :
  H = pyqubo.Constraint(
      sum(
          (X_i_t_m[i, t-1, m]  - X_i_t_m[i, 0, m] ) **2
          for m in K
          for i in CM),
          'w_1')
  return H

def build_rule8(X_i_t_m: pyqubo.array.Array) :
  H = pyqubo.Constraint(
      sum(
          (X_i_t_m[i, t-1, m]  - X_i_t_m[i, 0, m] ) **2
          for m in L
          for i in RM),
          'w_1')
  return H

def build_rule9(X_i_t_m: pyqubo.array.Array) :
  H = pyqubo.Constraint(
      sum(
          (X_i_t_m[i, t-1, m]  - X_i_t_m[i, 0, m] ) **2
          for m in Z
          for i in RM),
          'w_1')
  return H

def build_rule10(X_i_t_m: pyqubo.array.Array) :
  H = pyqubo.Constraint(
      sum(
          (sum(X_i_t_m[i, 0, m] for m in K) + sum(X_i_t_m[i, 0, m] for m in Z)) **2
          for i in GM),
          'w_1')
  return H

def build_rule11(X_i_t_m: pyqubo.array.Array) :
  H = pyqubo.Constraint(
      sum(
          (sum(X_i_t_m[i, t-1, m] for m in K) + sum(X_i_t_m[i, t-1, m] for m in Z)) **2
          for i in GM),
          'w_1')
  return H

def build_rule12(X_i_t_m: pyqubo.array.Array):
    H = pyqubo.Constraint(
        sum(
            (sum(X_i_t_m[i, t, m] for t in T for m in K) + sum(X_i_t_m[i, t, m] for t in T for m in Z) - 1) ** 2
            for i in GM),
        'w_1')
    return H

def build_rule13(X_i_t_m: pyqubo.array.Array) :
  H = pyqubo.Constraint(
      sum(
          (sum(X_i_t_m[i, 0, m] for m in L)) **2
          for i in CM),
          'w_1')
  return H

def build_rule14(X_i_t_m: pyqubo.array.Array) :
  H = pyqubo.Constraint(
      sum(
          (sum(X_i_t_m[i, t-1, m] for m in L)) **2
          for i in CM),
          'w_1')
  return H

def build_rule15(X_i_t_m: pyqubo.array.Array,X_i:pyqubo.array.Array) :
  H = pyqubo.Constraint(
      sum(
          (sum(X_i_t_m[i, t, m] for t in T for m in L) - X_i[i-g_num]) **2
          for i in CM),
          'w_1')
  return H

def build_rule16(X_i_t_m: pyqubo.array.Array,V_j:pyqubo.array.Array,d4:pyqubo.array.Array) :
  H = pyqubo.Constraint(
      sum(
          (sum(X_i_t_m[i, t, m]  for i in CGM) - V_j[m] + d4[0,m,t]) **2
          for t in T
          for m in K
      ),
          'w_1')
  return H

def build_rule17(X_i_t_m: pyqubo.array.Array,V_j:pyqubo.array.Array,e5:pyqubo.array.Array) :
  H = pyqubo.Constraint(
      sum(
          (sum(X_i_t_m[i, t, m]  for i in CRM) - V_j[m] + e5[0,m-k_c,t]) **2
          for t in T
          for m in L
      ),
          'w_1')
  return H

def build_rule18(X_i_t_m: pyqubo.array.Array,V_j:pyqubo.array.Array,f6:pyqubo.array.Array) :
  H = pyqubo.Constraint(
      sum(
          (sum(X_i_t_m[i, t, m]  for i in GRM) - V_j[m] + f6[0,m-k_c-l_r,t]) **2
          for t in T
          for m in Z
      ),
          'w_1')
  return H

def build_rule19(X_i_t_m: pyqubo.array.Array,V_j:pyqubo.array.Array,t) :
  H = pyqubo.Constraint(
      sum(
          (sum(X_i_t_m[i, t, m]  for i in CGM for t in T) - V_j[m]*(t+1)) **2
          for m in K
      ),
          'w_1')
  return H

def build_rule20(X_i_t_m: pyqubo.array.Array,V_j:pyqubo.array.Array,t) :
  H = pyqubo.Constraint(
      sum(
          (sum(X_i_t_m[i, t, m]  for i in CRM for t in T) - V_j[m]*(t+1)) **2
          for m in L
      ),
          'w_1')
  return H

def build_rule21(X_i_t_m: pyqubo.array.Array,V_j:pyqubo.array.Array,t) :
  H = pyqubo.Constraint(
      sum(
          (sum(X_i_t_m[i, t, m]  for i in GRM for t in T) - V_j[m]*(t+1)) **2
          for m in Z
      ),
          'w_1')
  return H
H = (build_objective(X_i,X_i_t_m,V_j)
     +\
      pyqubo.Placeholder('w_1') * build_rule1(X_i_t_m,X_i,a1)
     +\
    pyqubo.Placeholder('w_1') * build_rule2(X_i_t_m,X_i,b2)
     +\
    pyqubo.Placeholder('w_1') * build_rule3(X_i_t_m,X_i,c3)
     +\
    pyqubo.Placeholder('w_1') * build_rule4(X_i_t_m, V_j)
     +\
    pyqubo.Placeholder('w_1') * build_rule5(X_i_t_m,V_j)
     +\
    pyqubo.Placeholder('w_1') * build_rule6(X_i_t_m,V_j)
     +\
    pyqubo.Placeholder('w_1') * build_rule7(X_i_t_m)
     +\
    pyqubo.Placeholder('w_1') * build_rule8(X_i_t_m)
     +\
    pyqubo.Placeholder('w_1') * build_rule9(X_i_t_m)
     +\
    pyqubo.Placeholder('w_1') * build_rule10(X_i_t_m)
     +\
    pyqubo.Placeholder('w_1') * build_rule11(X_i_t_m)
     +\
    pyqubo.Placeholder('w_1') * build_rule12(X_i_t_m)
     +\
    pyqubo.Placeholder('w_1') * build_rule13(X_i_t_m)
     +\
    pyqubo.Placeholder('w_1') * build_rule14(X_i_t_m)
     +\
    pyqubo.Placeholder('w_1') * build_rule15(X_i_t_m,V_j)
     +\
    pyqubo.Placeholder('w_1') * build_rule16(X_i_t_m,V_j,d4)
     +\
    pyqubo.Placeholder('w_1') * build_rule17(X_i_t_m,V_j,e5)
     +\
    pyqubo.Placeholder('w_1') * build_rule18(X_i_t_m,V_j,f6)
     +\
    pyqubo.Placeholder('w_1') * build_rule19(X_i_t_m,V_j,5)
     +\
    pyqubo.Placeholder('w_1') * build_rule20(X_i_t_m,V_j,5)
     +\
    pyqubo.Placeholder('w_1') * build_rule21(X_i_t_m,V_j,5)
     )

feed_dict = {'w_1':5000,

             }
model = H.compile()
qubo, constant = model.to_qubo(feed_dict=feed_dict)

# 去掉 trotter 参数，或者替换为正确的参数名
sampler = SQASampler()  # 不传递参数，查看是否有其他方法设置 num_reads

response = sampler.sample_qubo(qubo)


# 假设 response 是采样器返回的结果
# 获取最优解
best_solution = response.first.sample  # 最优解的变量赋值
best_energy = response.first.energy    # 最优解对应的能量值

print("最优解：", best_solution)
print("最优解的能量值：", best_energy)


def evaluate_objective(X_i, X_i_t_m, V_j):
    H = (
        # 对于 CM 中的每个 i 计算
            sum(fix_c * X_i.get(f'X_i[{i - len(GM)}]', 0) for i in CM)
            + sum(fix_k * V_j.get(f'V_j[{i}]', 0) for i in K)
            + sum(
        C[i, j] * X_i_t_m.get(f'X_i_t_m[{i}][{t}][{m}]', 0) * X_i_t_m.get(f'X_i_t_m[{i}][{t+1}][{m}]', 0) * T_k
        for i in CGM for j in CGM for t in T1 for m in K if j != i
    )
            + sum(fix_r * X_i.get(f'X_i[{r - len(GM)}]', 0) for r in RM)
            + sum(fix_l * V_j.get(f'V_j[{m}]', 0) for m in L)
            + sum(fix_z * V_j.get(f'V_j[{m}]', 0) for m in Z)
            + sum(
        C[i, j] * X_i_t_m.get(f'X_i_t_m[{i}][{t}][{m}]', 0) * X_i_t_m.get(f'X_i_t_m[{i}][{t+1}][{m}]', 0) * T_l
        for i in CRM for j in CRM for t in T1 for m in L if j != i
    )
            + sum(
        C[i, j] * X_i_t_m.get(f'X_i_t_m[{i}][{t}][{m}]', 0) * X_i_t_m.get(f'X_i_t_m[{i}][{t+1}][{m}]', 0) * T_z
        for i in GRM for j in GRM for t in T1 for m in Z if j != i
    )
            + sum(
        C[i, j] * X_i_t_m.get(f'X_i_t_m[{i}][{t}][{m}]', 0) * X_i_t_m.get(f'X_i_t_m[{i}][{t+1}][{m}]', 0) * E_k
        for i in CGM for j in CGM for t in T1 for m in K if j != i
    )
            + sum(
        C[i, j] * X_i_t_m.get(f'X_i_t_m[{i}][{t}][{m}]', 0) * X_i_t_m.get(f'X_i_t_m[{i}][{t+1}][{m}]', 0) * E_l
        for i in CRM for j in CRM for t in T1 for m in L if j != i
    )
            + sum(
        C[i, j] * X_i_t_m.get(f'X_i_t_m[{i}][{t}][{m}]', 0) * X_i_t_m.get(f'X_i_t_m[{i}][{t+1}][{m}]', 0) * E_z
        for i in GRM for j in GRM for t in T1 for m in Z if j != i
    )
    )

    return H


# 假设解已经存储在字典中
# 从解中提取V_j和X_i_t_m的值
V_j = {key: value for key, value in best_solution.items() if key.startswith('V_j')}

X_i = {key: value for key, value in best_solution.items() if key.startswith('X_i[')}

X_i_t_m = {key: value for key, value in best_solution.items() if key.startswith('X_i_t_m')}

a1 = {key: value for key, value in best_solution.items() if key.startswith('a1')}

b2 = {key: value for key, value in best_solution.items() if key.startswith('b2')}

c3 = {key: value for key, value in best_solution.items() if key.startswith('c3')}

d4 = {key: value for key, value in best_solution.items() if key.startswith('d4')}

e5 = {key: value for key, value in best_solution.items() if key.startswith('e5')}

f6 = {key: value for key, value in best_solution.items() if key.startswith('f6')}
result = evaluate_objective(X_i,X_i_t_m,V_j)
print(result)

def evaluate_rule1(X_i_t_m, X_i, a1):
    H = (
        sum(
            (sum(X_i_t_m.get(f'X_i_t_m[{i}][{0}][{m}]', 0) for m in K) - lenk_c * X_i.get(f'X_i[{i - g_num}]', 0) + sum((2 ** l) * a1.get(f'a1[{l}][{i-g_num}]', 0) for l in range(a))) ** 2
            for i in CM))
    return H
result = evaluate_rule1(X_i_t_m, X_i, a1)
print(result)

def evaluate_rule2(X_i_t_m, X_i, b2):
    H =(
        sum(
            (sum(X_i_t_m.get(f'X_i_t_m[{i}][{0}][{m}]', 0) for m in L) - lenl_r * X_i.get(f'X_i[{i - g_num}]', 0) + sum((2 ** l) * b2.get(f'b2[{l}][{i-g_num-c_num}]', 0) for l in range(b))) ** 2
            for i in RM)
)
    return H
result = evaluate_rule2(X_i_t_m, X_i, b2)
print(result)

def evaluate_rule3(X_i_t_m, X_i, c3):
    H =(
        sum(
            (sum(X_i_t_m.get(f'X_i_t_m[{i}][{0}][{m}]', 0) for m in Z) - lenz_r * X_i.get(f'X_i[{i - g_num}]', 0) + sum(2 ** l * c3.get(f'c3[{l}][{i-g_num-c_num}]', 0) for l in range(c))) ** 2
            for i in RM))
    return H
result = evaluate_rule3(X_i_t_m, X_i,c3)
print(result)

def evaluate_rule4(X_i_t_m, V_j):
    H = (
        sum(
            (sum(X_i_t_m.get(f'X_i_t_m[{i}][{0}][{m}]', 0) for i in CM) - V_j.get(f'V_j[{m}]', 0)) ** 2
            for m in K))
    return H
result = evaluate_rule4(X_i_t_m, V_j)
print(result)

def evaluate_rule5(X_i_t_m, V_j):
    H = (
        sum(
            (sum(X_i_t_m.get(f'X_i_t_m[{i}][{0}][{m}]', 0) for i in RM) - V_j.get(f'V_j[{m}]', 0)) ** 2
            for m in L))
    return H
result = evaluate_rule5(X_i_t_m, V_j)
print(result)

def evaluate_rule6(X_i_t_m, V_j):
    H = (
        sum(
            (sum(X_i_t_m.get(f'X_i_t_m[{i}][{0}][{m}]', 0) for i in RM) - V_j.get(f'V_j[{m}]', 0)) ** 2
            for m in Z))
    return H
result = evaluate_rule6(X_i_t_m, V_j)
print(result)

def evaluate_rule7(X_i_t_m):
    H = (
        sum(
            (X_i_t_m.get(f'X_i_t_m[{i}][{t-1}][{m}]', 0) - X_i_t_m.get(f'X_i_t_m[{i}][{0}][{m}]', 0)) ** 2
            for m in K
            for i in CM))
    return H
result = evaluate_rule7(X_i_t_m)
print(result)

def evaluate_rule8(X_i_t_m):
    H = (
        sum(
            (X_i_t_m.get(f'X_i_t_m[{i}][{t-1}][{m}]', 0) - X_i_t_m.get(f'X_i_t_m[{i}][{0}][{m}]', 0)) ** 2
            for m in L
            for i in RM))
    return H
result = evaluate_rule8(X_i_t_m)
print(result)

def evaluate_rule9(X_i_t_m):
    H =(
        sum(
            (X_i_t_m.get(f'X_i_t_m[{i}][{t-1}][{m}]', 0) - X_i_t_m.get(f'X_i_t_m[{i}][{0}][{m}]', 0)) ** 2
            for m in Z
            for i in RM))
    return H
result = evaluate_rule9(X_i_t_m)
print(result)

def evaluate_rule10(X_i_t_m):
    H = (
        sum(
            (sum(X_i_t_m.get(f'X_i_t_m[{i}][{0}][{m}]', 0) for m in K) + sum(X_i_t_m.get(f'X_i_t_m[{i}][{0}][{m}]', 0) for m in Z)) ** 2
            for i in GM))
    return H
result = evaluate_rule10(X_i_t_m)
print(result)

def  evaluate_rule11(X_i_t_m):
    H = (
        sum(
            (sum(X_i_t_m.get(f'X_i_t_m[{i}][{t-1}][{m}]', 0) for m in K) + sum(X_i_t_m.get(f'X_i_t_m[{i}][{t-1}][{m}]', 0) for m in Z)) ** 2
            for i in GM))
    return H
result = evaluate_rule11(X_i_t_m)
print(result)

def evaluate_rule12(X_i_t_m):
    H = (
        sum(
            (sum(X_i_t_m.get(f'X_i_t_m[{i}][{t}][{m}]', 0) for t in T for m in K) +
             sum(X_i_t_m.get(f'X_i_t_m[{i}][{t}][{m}]', 0) for t in T for m in Z) - 1) ** 2
            for i in GM))
    return H
result = evaluate_rule12(X_i_t_m)
print(result)
def evaluate_rule13(X_i_t_m) :
  H = (
      sum(
          (sum(X_i_t_m.get(f'X_i_t_m[{i}][{0}][{m}]', 0) for m in L)) **2
          for i in CM))
  return H
result = evaluate_rule13(X_i_t_m)
print(result)
def evaluate_rule14(X_i_t_m) :
  H = (
      sum(
          (sum(X_i_t_m.get(f'X_i_t_m[{i}][{t-1}][{m}]', 0) for m in L)) **2
          for i in CM))
  return H
result = evaluate_rule14(X_i_t_m)
print(result)
def evaluate_rule15(X_i_t_m,X_i) :
  H = (
      sum(
          (sum(X_i_t_m.get(f'X_i_t_m[{i}][{t}][{m}]', 0) for t in T for m in L) - X_i.get(f'X_i[{i - g_num}]', 0)) **2
          for i in CM))
  return H
result = evaluate_rule15(X_i_t_m,X_i)
print(result)

def evaluate_rule16(X_i_t_m,V_j,d4) :
  H = (
      sum(
          (sum(X_i_t_m.get(f'X_i_t_m[{i}][{t}][{m}]', 0) for i in CGM) - V_j.get(f'V_j[{m}]', 0) + d4.get(f'd4[{0}][{m}][{t}]', 0)) **2
          for t in T
          for m in K
      ))
  return H
result = evaluate_rule16(X_i_t_m,V_j,d4)
print(result)
def evaluate_rule17(X_i_t_m,V_j,e5) :
  H = (
      sum(
          (sum(X_i_t_m.get(f'X_i_t_m[{i}][{t}][{m}]', 0)  for i in CRM) - V_j.get(f'V_j[{m}]', 0) + e5.get(f'e5[{0}][{m-k_c}][{t}]', 0)) **2
          for t in T
          for m in L
      ))
  return H
result = evaluate_rule17(X_i_t_m,V_j,e5)
print(result)
def evaluate_rule18(X_i_t_m,V_j,f6):
  H =(
      sum(
          (sum(X_i_t_m.get(f'X_i_t_m[{i}][{t}][{m}]', 0)  for i in GRM) - V_j.get(f'V_j[{m}]', 0) + f6.get(f'f6[{0}][{m-k_c-l_r}][{t}]', 0)) **2
          for t in T
          for m in Z
      ))
  return H
result = evaluate_rule18(X_i_t_m,V_j,f6)
print(result)

